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ANOTHER  STRATEGY  FOR  FAST  "POISSON" 

SOLVING  WITH  NON-CONSTANT  COEFFICIENTS 

1.  The  Introduction 

The  key  Co  fast  Poisson  solving  in  an  extended  domain  Is  nonlocal 

residual  error  dispersal.  On  a  two-dimensional  domain  characterized  by 

N  grid  points  In  one  direction  (taken  as  x  below)  and  M  grid  points  in 

another  (taken  as  y  ) ,  the  general  finite  difference  equation  could  be 
3 

solved  in  ~  (NM)  operations  by  matrix  inversion.  Because  the  matrix  is 
so  sparse,  however,  faster  algorithms  are  sought.  For  a  number  of  useful 
cases,  usually  characterized  by  cyclic  reduction,  Fourier  transforms,  or 
similar  folding  techniques,  fast  NM  log2NM  algorithms  exist  and  have  been 
in  use  for  some  time.  These  direct  solution  methods  are  non-iterative  and 
disperse  the  errors  globally  and  correctly  in  one  application  because  Che 
complete  elliptic  operator  can  be  inverted-. 

When  complex  stretched  grids  are  used  or  diffusion/dielectric  coef¬ 
ficients  vary  generally  in  space,  the  nature  of  the  physical  solutions  does 
not  change  but  the  variable  nature  of  the  finite-difference  coefficients 
in  the  resulting  elliptic  equations  obviates  the  use  of  direct  methods. 
Although  local  Inversion  approaches  such  as  ICCG,  SOR,  etc.  are  often  useful, 
there  are  still  good  reasons  to  seek  computationally  inexpensive  nonlocal 
techniques  whose  worst  case  convergence  is  expected  to  be  far  better  than  for 
local  inverse  techniques. 

A  nonlocal  technique  spreads  residual  error  at  each  point  all  over  the 

mesh  in  a  single  Iteration  cycle.  The  analytic  formulation  for  this  involves 
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the  manipulation  of  the  Green’s  function.  Numerically,  the  use  of  Green's 
function  in  a  general  case  can  be  prohibitively  expensive,  so  cheaply  manipu¬ 
lated  approximate  Green's  functions  are  sought.  It  is  not  necessary  for  an 
algorithm  to  be  local,  to  be  fast. 


II.  The  Idea 

This  short  note  contains  an  idea  for  fast  cheap  "Poisson"  solving  on 
a  rectangular  domain.  Consider  a  finite-difference  formulation  of  the 
equation 

V-7(x,y)  VP  (x,y)  -  »(x,y)  P(x,y)  -  S(x,y)  (1) 


vhere  v(x,y)  >  0,  u>(x,y)  a  Oand  (x,y)  in  the  domain(0  £  x  £  x  Osysy  _)  *• 

max,  *  'max 

Heat  transport,  special  treatments  of  fluid  dynamics,  and  electromagnetic 
propagation  in  varying  media  all  satisfy  similar  equations.  For  the  scalar 
field  |  (1  S  is  N,  1  s  j  S  H) ,  the  following  finite  difference  equation 

plus  well  posed  boundary  conditions  approximates  (1) , 

vi+l/2j<Pi+lj  "  Pij)  "  vi-l/2j  (Pij  “  pi-lJJ“  wi j  Pij  + 

y  y  (2) 

vij+l/2(Pij+l~  Pij*  "  vij-l/2(Pij  "  ^J-l*  "  Sij  * 

Suppose  is  unknown  with  {u^}  ,  jS^  }  ,  {  vj+1/2j  \  ,  and  8iven* 

A  new  set  of  unknowns  can  b*  defined  by  the  substitution  in  Eq.  (2), 

PiJ  "  *i/ij*  (i  "  1»  i  m  1 . M>*  C3) 


The  purpose  of  the  transformation  is  to  provide  a  set  of  adjustable  coefficients 
which  will  be  optimized  with  reference  to  the  coefficients|v*+j^2ji  * 

y 

ivij+l/2}  ,  |  iDjy }  but  not  with  respect  to  { S^|  which  may  change  from  one 
times tep  to  the  next. 
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By  appropiate  choice  of  *  £^e  coeff*c*-enta  ia  the  derived  Poisson- 

like  elliptic  equation  for{P^|  can  be  made  as  nearly  constant  as  possible. 
Since  constant  coefficient  problems  can  be  solved  using  fast  algorithms,  the 
non-constant  portion  of  the  coefficients  are  minimized  row  by  row  and  added 
as  a  residual  to  the  source  term  on  the  right  hand  side.  The  object  equation 
derived  from  (2)  by  substituting  (3)  is 

X  -X  y 

vi+l/2j  *i+lj  0i+lj  +  vi-l/2j  *i-lj0i-lj  +vij+l/2  *ij+l  *ij+l  + 
vij-l/2  *ij-llj-l  “  *ij(vi+l/2j  +  vi-l/2j  +  vij+l/2  +Vij-l/2  +  (4) 

"ip  0ij  "  Sij  * 

For  each  row  j  five  constants  (a^ ,  b^,  c^,  d^,  e^)  are  to  be  chosen  in 
conjunction  with  the  as  yet  undetermined  j  J  to  minimize  the  following  error 


measure: 


E  "j?i  “  “S  £ltaJ  ‘  *i+13l2  +  y?i  £lCl~  W 

jE  £itdj+i  ”  v^+i/2  +  €j?i£i  ej“i  ~  v^_i/2 


+s  s  ib.  -  ( — )±,  *±1r 

j-li-1  3  13  3 


where  I  have  designated  ( - ^  m  +  v*-i/2j  +  Vij+l/2  + 


vij-l/2  +u’ij)- 


The  undetermined  row  constants  ja^  | |  c j  |  >|d^  J  »  an<*  |  e j  }  '  ®l°ng  with  the 
undetermined  potential  multipliers  satisfy  a  set  of  equations  derived 


one  finds 


*ij  "  Fij  la  vi-l/2jaj+YVi+l/2jCj+6Vij-l/2dj+6vii+l/2ej+(  (8) 

where  Che  freedom  Co  choose  Che  oulciplicatlve  scale  of  in  Eq.  (3) 

pennies  us  Co  set  |  b^  |  -  1  without  loss  of  generality.  Here  I  have 

defined 


tu  .  +  y<»;+1/2/  +  ‘(v^-i/i)2  +  +  <— 1 '«>•  (9) 


The  row  constants  are  also  found  by  minimizing  E  .  Setting 
-  0,  gE/Sb  »  0,  etc.  gives  the  additional  four  equations 

i  i 

for  jej},  and  je^. 


Naj  “  S  v*+i/2j  Fi+lj  lQrVi+l/2j  aj  +  VVi+3/2j  Cj  + 

C10a) 

^1+13-1/2  dj  +  €vi+lj+l/2  *j  +  (  )i+lj1’ 

Ncj  Vi-l/2j  Fi-lj  Iavi-3/2  aj  +  Yvi-1/Z  CJ  + 

(10c) 

*i-lj-l/2  dj  +  €vi-lj+l/2  ®j  +  (  )i-lj1* 


D  X 

Nd  -S  Fij  lavi-l/2j  +  YVi+l/2*  CJ  + 

J  i«i 

**Ll/2  dJ  +  Hu+lrt  *1  +  (—V*  “d 

"•j  *  2^  »i1+i/2  Fu  >""1-1/21  *1  +  ''vl+l/2j  C1  + 

*"11-1/2  dl  +  €vW  *1  +  ’ ’«’• 


(lOd) 


(lOe) 
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These  four  linear  Inhomogeneous  equations  for  a^,  c^,  and  e^  can  be 
solved  without  reference  to  the  results  derived  from  any  other  row.  A 


block  tridiagonal  system  does  not  result. 


Given  these  row  coefficients,  which  can  be  seen  not  to  depend  on  the 
source  term  {s^} ,  Eq.  (8)yields  the  multiplicative  coefficients  { $  [  for 

each  of  the  unknown  potential  values  The  finite-difference  equation 


which  results  to  describe^  has  nearly  constant  coefficients,  row  by  row, 
the  complicating  spatial  variation  being  largely  cancelled  out  by  the 
derived  variations  of  }  . 

Ill .  The  Implementation 

This  algorithm  has  not  yet  been  implemented  numerically.  Using  the  row 

constants  derived  from  Eqs.  (10)  and  the  multiplication  coef f icientst^j } , 

the  left  hand  side  of  Eq.  (1)  can  be  expanded  as  follows: 

£•  v(x,y)  7  l^x.y)  -  ui(x,y)  P(x,y)  -  D2  0(x,y)  +  [7'v  7  7  *  (*,y) 

-  m  1  *  (x,y)  -  D2  0  (x.y)  ).  (11) 

2 

The  operator  D  ,  by  construction,  is  Fourier- transformable  in  the  x 


direction  and  results  in  decoupled  tridiagonal  equations  in  the  y  direction. 

2 

The  second  term  on  the  right  in  (11),  which  I  shall  call  R  ,  has  been 
minimized  in  a  simple  least  squares  sense.  Rewriting  Eq.  (1)  gives 

"  (x.y)  0 (x.y)  - *2(x,y)  -  R*(x,y)  *®(x,y)  (12) 

as  an  iterative  form  for  solving  (1)  which  should  converge  rapidly.  The 
superscript  a  indicates  which  level  of  Iteration  has  been  passed.  The 
right  hand  side  of  (12)  is  initialized  using  0°(x,y),  a  guess  at  the 
expected  potential. 
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